ASYMPTOTIC-NUMERICAL STUDY OF SUPERSENSITIVITY FOR 
GENERALIZED BURGERS EQUATIONS 
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Abstract. This article addresses some asymptotic and numerical issues related to the solution of 
Burgers' equation, —eu xx -\-ut+uu x = on (—1, 1), subject to the boundary conditions u(— 1) = 1+5, 
= —1, and its generalization to two dimensions, — eAu+ut + uu x + uu y = on (—1, 1) X (— n, it), 
subject to the boundary conditions u\ x= i =1 + 5, u\ x= -i = —1, with 2n periodicity in y. The 
perturbation parameters 5 and e are arbitrarily small positive and independent; when they approach 
0, they satisfy the asymptotic order relation 5 = 3 (e~ a / e ) for some constant a £ (0, 1). 

The solutions of these convection-dominated viscous conservation laws exhibit a transition layer 
in the interior of the domain, whose position as t — > oo is supersensitive to the boundary perturbation. 
Algorithms are presented for the computation of the position of the transition layer at steady state. 
The algorithms generalize to viscous conservation laws with a convex nonlinearity and are scalable 
in a parallel computing environment. 
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1. Introduction. In this article we address some asymptotic and numerical is- 
sues related to the solution of Burgers' equation, 

(1) - eu xx + u t + uu x = on (-1,1), u(-l) = l + <5, u(l) = -1, 
and its generalization to two dimensions, 

(tyeAu + ut + uu x + f3uu y = on (— 1, 1) X (— tt, tt), u\ x =~i = l + <5, u\ x =i = — 1. 

In the latter case, we assume periodicity (period 2ir) in the second coordinate (y). The 
perturbation parameters S and e are arbitrarily small positive; they arc independent, 
but when they approach 0, they satisfy the asymptotic order relation 

(3) 5 = O s {e~ a/£ ) as<5,e|0, 

for some constant a £ (0, 1) which does not depend on S or e. This asymptotic relation 
implies that 5 is transcendentally small (in the sense of asymptotic analysis) compared 
with e, but S dominates e -1 / £ as e i 0. (Sec g [| |, f§ § for definitions and basic 
concepts of asymptotic analysis.) 

If £ = 0, the solution of ([!]) develops a shock (discontinuity) in finite time, even 
when the initial data are smooth || The perturbation introduced by a nonzero 
e models the presence of viscosity, which tends to smooth the discontinuity ||, |J. 
Instead of a shock, one has a transition layer — a region of rapid variation, which 
extends over a distance 0(e) as e [ 0. The position of the transition layer varies 
with time, and its eventual location at steady state is extremely sensitive to the 
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boundary data. In fact, even the transcendentally small perturbation S leads to a 
measurable (that is, order one) effect on the eventual location of the transition layer. 
This phenomenon, known as supersensitivity, was first observed by Lorentz It 
has been studied extensively for Burgers' equation and more general viscous conser- 
vation laws in one dimension by Kreiss and Kreiss |]lT|, Kreiss fl2|j , Laforgue and 
O'Malley || @, || |f , 0, ||, and Reyna and Ward[(L9[ |cj |l| . 

An example from combustion theory shows that supersensitivity is of more than 
mathematical significance. A simple model of flame propagation in gaseous fuels 
involves a system of two coupled convection-diffusion equations, one for the tempera- 
ture of the mixture, another for the concentration of the reaction-limiting component 
in the mixture ]22| , § 3.2]. If one ignores exponentially small perturbations in the 
data, one finds that the Lewis number C, which is a measure for the ratio of heat 
and mass transfer in the mixture, has no effect on the location of the combustion 
front. Yet, numerical computations show that this location is very sensitive — in fact, 
supersensitive — to the value of C. 

Although the phenomenon of supersensitivity is fairly well understood theoreti- 
cally, at least for one-dimensional problems, the numerical solution of such problems 
still poses formidable challenges, especially in more than one dimension. The meth- 
ods that have been proposed in the numerical literature for singularly perturbed 



boundary- value problems (see, for example, |23 ) tend to focus on uniform approxi- 
mations or finite elements with special features, not on the supersensitive dependence 
of the transition layer on the boundary data. On the other hand, the algorithms 
we propose are designed specifically to capture the phenomenon of supersensitivity. 
They use the fact that the solution approaches a certain profile as the perturbation 
parameters approach zero and focus on the computation of the corrections. 

Our ultimate goal is to develop algorithms for multidimensional problems that 
are, first of all, suitable for long-time integration, so stable steady states can be com- 
puted with confidence; second, extremely accurate in space, so the eventual location 
of transition layers can be predicted with accuracy; and third, scalable in a multi- 
processing environment, so large-scale problems can be solved in a reasonable length 
of time. Although we discuss only Burgers' equation and its generalization to two 
dimensions, the algorithms are not restricted by the special form of the nonlincarity. 

In § 2 we consider Burgers' equation. We propose a simple algorithm that effec- 
tively captures the supersensitive location of the transition layer at steady state. We 
stress the importance of the regions outside the transition layer, where the solution 
does not yet deviate appreciably from the boundary values. In § 3 we address the 
generalized Burgers equation in two dimensions. We show through a formal asymp- 
totic analysis that the location of the transition layer may vary in the direction of 
periodicity (y), but the transition layer is essentially flat, and only its average position 
(averaged over y) depends supersensitively on the small parameters. We then develop 
an algorithm that effectively approximates the transition layer. 

2. One-Dimensional Case. We begin by considering Eq. (|l|), 

(4) - eu xx + u t + uu x = on (-1,1), u(-l) = 1 + S, u(l) = — 1. 

As shown by Laforgue and O'Malley jig , the solution approaches a certain profile 
function as e J, 0, 

(5) u(x, t) = tanh?? + e _ °/ e Mi(ry, a) + . . . , 
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where 

(6) r]=X _x^ 1 ^ a = te -«/c. 

e 

The hyperbolic tangent incorporates a transition layer centered at x* , which connects 
the limiting values ±1 at =Foo. Note that these limiting values are transcendentally 
close to the prescribed boundary values 1 + 6 and —1 of u at —1 and 1. The position 
of the transition layer varies on a transcendentally slow time scale; if <5 = 26e~°/ e , its 
limit as a — > oo is 

(7) < s = l-a + eln6. 

The important points to observe are that, first, the solution u approaches a certain pro- 
file function as e j 0; second, the accurate determination of the steady-state position 
of the transition layer requires long-term integration; and, third, a transcendentally 
small perturbation of the boundary data has a measurable effect on the location of 
the transition layer. 

The asymptotic analysis has been generalized to more general nonlinearities by 
Laforgue and O'Malley |l7], [l8| and Reyna and Ward Q, with similar results. One 
finds a limiting profile, which generalizes the hyperbolic tangent function, and a transi- 
tion layer which moves on a transcendentally slow time scale to a steady-state position. 
This position depends supersensitively on the boundary perturbation. Whenever this 
situation arises, appropriate variants of the following algorithms can be developed. 

2.1. Spatial Approximation. To approximate the solution in space, we use 
a domain decomposition method with two non-overlapping subdomains, where the 
interface is located approximately at the center of the transition layer, an adaptive 
pseudo-spectral method on each subdomain, and collocation based on Tchebychev 
polynomials, where the collocation points are concentrated in the transition layer. 
The algorithm is standard and has been described elsewhere [g4|, |25|, ^7], |2q |; we 
summarize it here only for completeness. 

Let x* S (— 1, 1) denote the (approximate) position of the center of the transition 
layer; x* varies in time (t), but since t enters only as a parameter in the discussion of 
the spatial approximation, we do not write it explicitly. We decompose, 

(8) ^ = (-1,0, n 2 = (a;*,l) ) 

and map each of the subdomains Qi and linearly onto (—1, 1), 

gx : y e (-1, 1) h- x = gi(y) = -1 + §(z* + l)(y + 1) E fi x , 

g 2 :ye (-1, 1) h- x = g 2 (y) = 1 - i(l - ar*)(l - y) S il 2 . 

The restrictions of u to f^i and exhibit boundary layer behavior near x* . The point 
x = x* corresponds to y = 1 under g\ and to y = —1 under 172- To concentrate the 
collocation points near x*, we define a one-parameter family of nonlinear mappings 
of the interval (—1, 1) onto itself, 

/i(- ,a) : s E (-1, 1) y = /i(s,a) = 1 - (4/7r)arctan (atan \(l - s)tt) e (-1, 1), 



h(- ,a) : s e (-1, 1) i-v y = / 2 (s, a) 



1 ■ (4/7r)arctan (atani(s + l)ir) <E (-1, 1). 
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If the parameter a is small, fx concentrates points near 1 and / 2 concentrates points 
near —1. Concentrating points near critical points is the computational analog of 
coordinate stretching in asymptotic analysis. The choice of a can be optimized by 
means of a priori estimates [^J, p5| ; we usually take a = e 1 ! 2 p7| . The composite 
maps, 

hi(-;a)=g i (-)of i (-,a), i = l,2, 

are one-to-one from (—1,1) onto Qf, we denote their inverses by ; a), i = 1, 2. 

We look for the solution of Eq. (|J) by approximating locally on each of the sub- 
domains fli and Q2 and imposing C 1 continuity at x*. If U denotes the global 
approximation, then 

U(x) = Ui{x) for a: £ O,, U G C 1 ([— 1, 1]). 

The local approximations consist of finite sums of Tchebychev polynomials, 

JV-l 

Ui(x) = ".//;'/', a)), x e Qi, i = 1,2; T^cosO) = cos(j0), 9 = ir/N. 
3=0 

2.2. Integration for Times of Order One. For short-time integration it suf- 
fices to combine a one-step forward Euler approximation with an implicit treatment 
of the second-order spatial derivative and an explicit treatment of the nonlinear term. 

Starting with an approximate solution U° = U{- ,<o) a t time toj w e identify the 
point x* = x*(to) with the location of the zero of U , partition the domain in two 
subdomains, and select the collocation points. Fixing this configuration temporarily, 
we compute a sequence {U n : n = 1, 2 . . . } of successive approximations U n using the 
algorithm 

77™ _ Tjn-i 

(9) - eD 2 U n + — + U n - x DU n - x =0, n = 1, 2, ... . 

The symbol D represents the pseudo-spectral differentiation operator in physical 
space |^9| . The time step At is constant, so U n is the approximate solution of the 
boundary-value problem (|l|) at £o + nAt. The algortithm (^|) is nonconservative. 

When the location of the zero of U n has moved over a distance e, we suspend the 
algorithm (|^). We shift x* to the current location of the zero, reconfigure the partition, 
update the collocation points, replace U° by the values at the new collocation points 
(using interpolation if necessary), and continue the algorithm (^. We repeat this 
process until the steady state is reached. The change in the location of the zero of the 
computed approximation becomes smaller as time progresses, so the same collocation 
configuration serves for longer and longer time intervals. 

The algorithm (^) requires the solution of U n from the equation 

(10) AU n = IT 1 - 1 + {A^U^DU"- 1 . 

The matrix A, which is order 27V — 1, has a block structure, 

I A x h \ 
A = \ a\ c b\ . 
V a 2 A 2 j 
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Ai and A 2 are square matrices of order N — 1; a\, a 2 , b\, and 6 2 vectors of length 
N—l; c is a constant. The center row accounts for the C 1 continuity at the interface. 
This block structure allows a solution of the system (|l^) in two parallel processes from 
opposite ends. The matrix A does not change as long as the collocation configuration 
is frozen. However, its condition number increases with the number of collocation 
points. This increase puts a lower limit on the values of S one can handle in practice. 



Table 1 

Location of the transition layer at steady state. 





£ = 0.1 


e = 0.05 


6 










1.0 ■ 10" 1 


0.72464 


0.700427 


0.86237 


0.850213 


1.0 ■ 1Q-' 2 


0.47486 


0.470176 


0.73755 


0.735084 


1.0 ■ lO" 3 


0.24133 


0.240724 


0.62055 


0.619955 


1.0 ■ 10~ 4 


0.05265 


0.052606 


0.50485 


0.504826 


1.0 • 10~ b 


0.00537 


0.005504 


0.38962 


0.389696 



Table 1 gives the location of the transition layer at steady state, x^, computed 
with the algorithm ([)]) with N = 39 collocation points in each subdomain and a time 
step At = 0.02. The number x* s = 1 — £ln(2/<5), which is an asymptotic estimate of 
x* (see Eq. (Q)) is given for comparison. The initial conditions were usually obtained 
by linear interpolation from the boundary data, but variations were made to test the 
answers. The lower limit on e is determined by the fact that the computation time for 
the algorithm (^) increases as e decreases. In § 2.4 we discuss an algorithm suitable 
for long-time integration. 

Table 2 shows the impact of grid refinement (the number of collocation points, 
N) on the value of x ^ . 



Table 2 

Effect of grid refinement (N ) on x*^; £ = 0.1, <5 = 1.0 ■ 10 -3 . 



N 


15 


19 


29 


39 


49 


59 




0.25576 


0.24115 


0.24166 


0.24133 


0.24140 


0.24143 



2.3. Neglecting Viscosity. In supersensitive boundary-value problems, the so- 
lution in the "tails" on either side of the transition layer is exponentially close to the 
prescribed boundary values, so the viscous term is exponentially small there. It is 
therefore tempting to assume that one can sacrifice some accuracy in the computa- 
tion of the viscous term during the transient phase and still find the position of the 
transition layer at steady state with a high degree of accuracy. 

An extreme form of this assumption underlies the approach where one constructs 
a first approximation by ignoring the viscous term altogether. Using a conservative 
finite-difference scheme, such as Godunov, one constructs the entropy solution of the 
inviscid conservation law (e = 0) and takes this as a first approximation. One then 
constructs higher-order uniform approximations, for example by means of a hetero- 
geneous domain-decomposition method, using either a different numerical scheme to 
solve the full viscid problem in the interior of the transition layer or some other 
approximation of the viscous equation. 

The x method introduced by Brezzi et al. |3C[| is a more sophisticated nonlinear 
adaptive scheme based on the same assumption. Here, one replaces the boundary- 
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value problem by 

(11) - e X (u xx ) + u t + {f(u)) x - on (-1, 1), u(-l, t) = l + 5, u(l, t) 



-1, 



where \ = X<j,t is a smooth monotone function, x(s) = if s < a and x( s ) = s 
if |s| > a + t for some positive numbers a and r. This method has been applied to 
Burgers' equation (3l) and the incompressible Navier-Stokes equations [j32| . However, 
we claim that the \ method cannot accurately predict the ultimate position of the 
transition layer, at least for Burgers' equation on a finite interval with Dirichlet data. 
This claim is supported by the following observations. 

Consider the results quoted in ]3l|, Table II]. With few exceptions, they involve 
relatively large values of a {a is called S in |H]]); in fact, a is typically greater than 
e -1 / 2 (e is called v in |3TJ ) by one or two orders of magnitude. The viscous term 
is therefore always neglected, unless u xx is of the same order as e -1 / 2 ; that is, the 
viscous term is neglected everywhere except in the transition layer. If the x method 
gave the correct position of the transition layer at steady state, then the same would 
certainly be the case when we simply multiply the viscous term by a smooth function 
of position, whose support is of order one and includes the transition layer. After all, 
in the latter case we account for the viscous term over a much broader region. These 
arguments lead us to consider the boundary-value problem 

(12) - eH(x)u xx + u t + uu x = on (-1, 1), tt(— 1, t) = 1 + 6, u(l, t) = -1, 

instead of the boundary- value problem (|ll"l). Here, H is a smooth cut-off function, 



(13) 



_ } h (1 - tanh(a(a; - x*(t) + /?))) , x<x*(t), 
W "l | (l-tanh(a(-z + +/?))), a? > ar*(t). 



We use a numerical approximation of x* (t) and choose the parameters a and j3 so 

H(x) = 1 if \x - x*(t)\ < H(x) = if \x - x*(i)\ > §/?, 

to within machine accuracy (1 • 1CT 15 ). The algorithm (|) leads to the solution of U r ' 
from the equation 



(14) - eHD 2 U n 



U n - U r ' 



U'^DvJ 1 - 1 = 0, n=l,2,. 



The results given in Table 3 (computed for e = 0.1 and S = 1.0 • 10~ 2 , with a = 200) 
show that the algorithm ( |l4| ) can give incorrect results for the position of the transition 
layer at steady state. (The correct value is x*^ — 0.47486, see Table 1.) The position 

Table 3 

Location of the transition layer at steady state computed with the \ method. 





0.1 


0.2 


0.3 


0.4 


0.5 


0.6 


0.7 




0.0056 


0.0110 


0.0174 


0.0384 


0.0923 


0.1067 


0.1565 



of the shock freezes too early during the transient phase. The actual moment of 
freezing depends on the size of the zone to the left of the transition layer (where 
H = 0). The result improves as (3 increases, but for f3 = 0.8 the algorithm fails to 
converge. (The position of the transition layer keeps oscillating between the regions 
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where H = 1 and H = 0.) On the other hand, if we include the missing part of the 
viscous term explicitly and use the algorithm 

77™ — /T™ -1 

{\h)eHD 2 U n - e(l - iJ)!? 2 ?/"- 1 + — + U^Du* 1 ' 1 = 0, n = 1,2, . . . , 

instead of (|l4|), we retrieve the correct position of the transition layer at steady state. 
This result demonstrates clearly that, when the problem is supersensitive, it is not 
advisable to neglect the viscous term, even when that term is exponentially small. 

Note that the modified algorithm ([D]) treats the second-order derivative explic- 
itly in the region where H = and implicitly in the region where H = 1. The idea 
of using a cutoff function to construct a composite algorithm is described in detail in 
our article [E8|. The procedure offers a very general tool for the design of heteroge- 
neous domain decompositions in the framework of a finite-difference approximation. 
However, since the algorithm ( |l5|) is based on a Tchebyshev pseudo-spectral approxi- 
mation, the partially explicit treatment of the viscous term forces a severe constraint 
on the time step that cannot be circumvented. For example, with e = 0.1 and N — 49 
collocation points per subdomain, the time step must be 10 times smaller than the 
time step for the algorithm (^). The algorithm (|l5|) is therefore certainly not practical 
for long-time integration. 

2.4. Long-Time Integration. The explicit treatment of the nonlinear term 
in the algorithm (Q) imposes a severe constraint on the time step (CFL condition). 
The algorithm is therefore not suitable for long-time integration. The alternative 
approach commonly taken is to use a fully implicit scheme in combination with a 
Newton algorithm j33|. However, an implicit scheme can be expensive and is certain 
to increase interprocessor communication in a multiprocessing environment. If the 
solution of the viscous conservation law is close to a certain profile function, as is the 
case for Burgers' equation, at least after an initial transient, the following algorithm 
offers a more efficient alternative. 

We start the integration of Eq. (^) at t = to, say, when we have an approximate 
profile with a transition layer centered at x* . We construct the function uq, 

x — x* 

(16) u (x) = — tanh 



2e ' 

which satifics Burgers' equation exactly, and look for a solution u of the form 

(17) u(x, t) = uo(x) + Sv(x, t). 
Then v must satisfy the nonlinear boundary- value problem 

-ev xx + v t + u v x + u' Q v + 5vv x = on (-1, 1), 

v(-l,t) =S- 1 (l + S-u Q (-l)), v(l,t) = 6- 1 (-l-u (l)). 
We integrate this boundary- value problem for t > t , using the algorithm 

yn _ yn-1 

(18) -£D 2 V n + — +u DV n + u' V n = -5V n - 1 DV n ~ 1 , n= 1,2,... , 

and define an approximation U of u for t > t , 



(19) 



U(x,t) = uo(x) + 5V(x,t), t > to- 
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We proceed with the integration as long as the supremum of U{- ,t) — u remains 
of the order of S. When this criterion is no longer met, at t = t\ say, we suspend 
the integration, identify the point x* with the location of the center of the transition 
layer at t\, update the profile function uq, and continue the integration beyond t\. 
We repeat the procedure until the steady state is reached. Figure 1 shows a profile 
function uq and the computed solution U(- , f) at some time t. Note that the former 
is monotone, the latter is not. 

1.5 

l 

0.5 

-0.5 

-1 
-1.5 

-2 

-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 

x 

Fig. 1. The profile function uo (dashed line) and the computed solution U(- ,t) at some time 
t (solid line); e = 0.05, 5=1- 10" 3 . 

The algorithm @ is very similar to (^|). The spatial approximation is handled 
with the same domain-decomposition method, so the structure of the resulting linear 
system is the same as in Eq. ( p^|). But the constraint on the time step is relaxed by 
a factor 5. With the algorithm (|18|) we can integrate the boundary-value problem for 
values of e as small as 0.01, down to Ss of the order of 10~ 6 , using time steps that 
are typically 50 times larger than with the algorithm ([^). Table 4 gives some results 
for a;^,,, computed with the algorithm (|l^), with N = 39 collocation points in each 
subdomain. 



Table 4 

Location of the transition layer at steady state. 





£ = 0.1 


e = 0.05 


e = 0.02 


£ = 0.01 


8 


x oc 








1.0 ■ 10" 1 


0.72346 


0.86175 






1.0 ■ 10^ 


0.47508 


0.73774 


0.89518 




1.0 • lO" 3 


0.24140 


0.62057 


0.84827 


0.92429 


1.0 ■ 10~ 4 


0.05275 


0.50485 


0.80210 


0.90104 


1.0 • 10" s 


0.00561 


0.38964 


0.75609 


0.87800 


1.0 ■ 10 _B 


0.00066 


0.27452 


0.70996 


0.85482 


1.0 ■ io-'' 








0.83084 
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3. Two-Dimensional Case. Next, we consider Burgers' equation generalized 
to two dimensions, 

(2Q$Au + u t + uu x +f3uu y = on (— 1, 1) X (— 7r, it), u\ x= -i = 1 + 5, u\ x= i = — 1. 

We assume periodicity in y (period 2tt). The perturbation 5 may vary with y; its 
Fourier expansion is 



The coefficients r/k, as well as the pre- factor So, are independent of y. The pre-factor 
Sq is arbitrarily small positive and defined in such a way that rjo — 1 and r/k — 0(1) 
as 5q I for k = ±1, ±2, .... The order relation (||) between S and e, which must 
hold uniformly in y, implies that 



for some a G (0,1). We show in § 3.1 that, under these conditions, the average 
position of the transition layer (averaged over y) depends supersensitively on the 
small parameters e and Sq- In § 3.2 we briefly review the spatial approximation. § 3.3 
is devoted to the integration procedure. 

3.1. Profile Function. The algorithm we propose for the solution of the boundary-^ 

value problem (|2^) is again based on the assumption that the solution is close to a 
known profile function. Our goal in this section is to show that, under the conditions 
given above, the profile function is asymptotically independent of y and again given 
to leading order by the hyperbolic tangent, as in Eq. (^). 

We introduce the constant x* € (0, 1) such that x* ~ 1 — eln(2/<5o) as £ J. 0. We 
define the function uq, 

x — x 

(23) uq{x) = — tanh — , 

and look for a profile function ip of the form 



(21) 




(22) 



So = O s (e- Q/e ) as<5 , £ j0, 



(24) 



p{x,y) = u (x) +v{x,y). 



Because uq satisfies Burgers' equation, v must satisfy the differential equation 



(25) 



eAv + uqv x + u' v + 0v,QVy + vv x + j3vv y = 0, 



together with the boundary conditions 



(26) 



S-S { 



v 



x=l 



= 0. 



The linearized equation, 



(27) 



£(v) 



eAv + uqv x + u v + fiuoVy = 0, 



has a solution, £(u' ) = 0, so we reduce the order by substituting 



(28) 



V = UqW. 
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Then w must satisfy the equation 

(29) — sAw — uqw x + (3uow y = 0, 
together with the boundary conditions 

(30) w\ x= -i = (5 - 5 )K(-1), H*=i = 0. 

We estimate w by means of the coefficients in its Fourier expansion, 

(31) w(x,y) = J2^k(x)e iky . 

kez 



The leading coefficient wo is the solution of the boundary- value problem 

(32) - ew'o - u w' = on (-1,1), w (-l) = 0, w (l) = 0. 

so wo = 0. Note that this result is a direct consequence of the fact that we have 
defined x* in terms of the average 5q in Eq. ((2^); any other definition leads to an 
inhomogeneous boundary-value problem, whose solution wq does not vanish. 

The remaining coefficients w k , k — ±1,±2,..., are found from the boundary- 
value problem 

(33) - sw k - u w' k + (ek 2 + i/3ku )w k — 0, Wk(-1) = SoVk, w k {l) = 0. 

This is a classical turning-point problem, as uq changes sign in the interval (—1,1). 
The asymptotic behavior of w k as e J, can be found by the method described in 
§3.E], 

{M)w k {x) - \c k e i0kx + (l - c k e~ lf)k ) e - {1+x}/e + (0 - c k e lf3k ) e - (1 ~ x)/e ] w k (-l). 
The coefficient c k is such that the functional 

(35) C[w] = [ew 12 + (ek 2 + i(3ku )w 2 ] exp j u (£,) dA dx 
has a critical point at w = w k , 

(36) ^1 = 0. 

dc k 

Notice that the differential equation (|3^) is the Euler equation of the functional £; 
the Neumann boundary conditions u>J..(±l) = are the natural boundary conditions 
associated with C. 

Obviously, finding an explicit expression for c k is out of the question. The best 
we can aim for is an asymptotic expansion as e j 0, and even here we must resort 
to computational assistance. Using the symbolic manipulation language MAPLE, we 
find 

(37) Cl ._f±^L e -.«.«.-, as6i0 . 

The first term in the brackets in Eq. ( |34|) represents the regular part of the asymptotic 
behavior of w k , which dominates in the interior; the remaining two terms represent 
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the singular part, which dominates near the endpoints of the interval. Since we are 
interested in the transition layer, which is located in the interior, we ignore the singular 
part and take only the regular part, Wk(x) ~ CkWk(— l^ 4 ' 3 '™. That is, we take 



SoVkC u' (x) pimi _ 2x * +x) 

(l + /3 2 )e 2 P u' (-l) 



(38) ^-TT^^^ 1 as ^°- 



If we use the asymptotic approximation 

u' (x) _ 1 -tanh 2 ((a;-a;*)/(2e)) 1 - tanh 2 ((x - x*)/{2e)) 
«o(-l) ~ l-tanh 2 ((-l-a;*)/(2£)) 4c -(i+^)/e ' 

we obtain the asymptotic expression 

This result implies that the Fourier series off, as well as those of v x and w y , converge. 
Furthermore, II are O(5 s 2 ). 

Finding the asymptotic behavior of ||ue||oo is less obvious. It follows from Eq. ( |3~4| ) 
that w' k (x) ~ Cke l P kx as e j 0, at least for cc in the interior of the domain. Therefore, 
v' k (x) = (u' w k Y(x) ~ (u'£w k )(x) ~ e _1 (u^w fe )(a;) = e -1 u fc (x). Hence, Halloo = 
O(6 £~ 3 ) as (5 ,£ | 0. 

Since S = O s (e- a / e ) for some a E (0, 1), we have S e~ p = O s {e~ a '' e ) (p = 2,3) 
for any a' 6 (0, a), so any solution w of Eq. ( p7| ) which satisfies the boundary con- 
ditions (26) is transcendentally small. The residue vv x + j3vv y , which was ignored 
in the transition from the nonlinear equation ( p5| ) to the linear equation ( p7| ) is like- 
wise transcendentally small and, in fact, 0(<5q£ -5 ), so we also have an a posteriori 
justification for the linearization. 

These arguments motivate the choice of uq, which depends only on x, as the 
profile function in the design of the numerical algorithms for Eq. (p0[). 

3.2. Spatial Approximation. The spatial approximation is again based on a 
domain decomposition with two non-overlapping subdomains on either side of the 
y-averaged location of the center of the transition layer. On each subdomain we use 
an adaptive pseudo-spectral method in the x direction and a finite-difference method 
in the y direction. The pseudo-spectral method is the same as in the one-dimensional 
case; it uses Tchebychev polynomial collocation with N x collocation points. 

Since the transition layer is close to a plane parallel to the x axis, there is no 
need to resort to an adaptive grid in the y direction. For our numerical experiments 
we chose a regular grid with mesh width h = 2ir/N y and a sixth-order central finite- 
difference approximation of u yy and u y . The choice may seem inconsistent with the 
spectral approximation in the x direction; a more obvious choice would be a Fourier 
approximation in the y direction. Theoretically, the finite-difference approximation 
in the y direction restricts the accuracy of the approximation for a regular problem 
(e = O s (l)) to sixth order, less than the accuracy guaranteed by the pseudo-spectral 
approximation in the x direction. However, as the transition layer is close to a plane 
parallel to the x axis, it is relatively easy to keep the numerical error in the finite- 
difference approximation of the term eu yy smaller than the numerical error in the 
pseudo-spectral approximation of the term eu xx with a moderate number of dis- 
cretization points N y . There is, therefore, no need to use a better approximation, like 
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the Fourier approximation, for the term su yy . Furthermore, the spectral radius of D 2 
is smaller with sixth-order finite differences than with Fourier differentiation. This 
difference implies an additional advantage for a finite-difference approximation when 
the y derivatives are treated explicitly |3jJ] . 

3.3. Integration for Times of Order One. We extend the Euler scheme @ 
to two dimensions as follows: 
Tjn _ 771—1 

-eD 2 x U n + — =eD 2 y U n - 1 + U n - 1 D x U n - 1 +pU n - 1 D y U n -\ n= 1,2,... . 

(40) 

Here, D x is the pseudo-spectral differential operator with Tchebychev polynomials, 
D y the finite-difference operator with sixth-order central finite differences. 

The algorithm ([13) has several features that make it readily parallelizable. First, 
the approximations Dyll 12 " 1 and DyU n ~ x are taken explicitly, so the variable y is 
only a parameter. Second, because we are using finite-difference approximations, we 
have only local data dependencies. This latter point especially offers a significant 
advantage over a spectral method, which uses global interpolation. 

Table 5 shows the results for the boundary- value problem (|20| ) with j3 = 1, e = 0.1 
and piecewise constant boundary data with So = 1.0 • 10~ 2 , 

( 1.01 -AS if -7r < y < -^71-, 
(41) u(-l,y) = < 1.01 + AS if - ±tt < y < ±tt, 

[ 1.01 -AS if ±7T < y < 7T. 

The algorithm ( |40| ) was applied with = 39 collocation points per subdomain in 
the x direction and N y = 32 interpolation points in the y direction. The table gives 
the y-averaged location of the transition layer at steady state, <x*„ >, as well as the 
maximum deviation of the center of the transition layer from its y-averaged location, 
Ax*; that is, x* yo (y) varies between <x* c > — Ax* and <x* XJ > + Ax*. These results 

Table 5 

Location of the transition layer at steady state; boundary data (U4). 



A<5 


<X*oo> 


Ax* 


0.25 • 10~* 


0.4758 


1.3114 • 10-* 


0.50 ■ 10-* 


0.4759 


2.3584 ■ 10-* 


1.0 ■ 10 _a 


0.4758 


4.8865 • 10-* 


1.5 ■ 10"^ 


0.4750 


7.8632 • 10~* 


2.0 ■ 10"^ 


0.4738 


9.3250 ■ 10~* 


3.0 ■ 10-* 


0.4700 


15.688 • 10 _:i 



show that the algorithm ( |40| ) is extremely effective for the boundary-value problem. 
As AS increases, Ax* grows approximately linearly with AS. The graph of U — Uq 
maintains its overall shape, so its width (which measures Ax*) varies in proportion 
to its height (which measures \\U — uo||oo). The numerical results therefore indicate 
also that ||{7 — ito[[oo grows approximately linearly with AS. This conclusion matches 
the results of the asymptotic analysis in § 3.1, in particular Eq. (|39|), where it was 
shown that Vk is proportional to 77^. 

Results for a much harder case are presented in Table 6. The parameters (3 and 
e are fixed as before, (3 — 1 and e = 0.1, but this time the data at the left boundary 
are sharply peaked at the midpoint, 

(42) u(-l,y) = 1.005+ (A^e- 20 ^-™ 8 ^, -tt < y < n. 
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The algorithm ( [40| ) was again applied with N x — 39 collocation points per subdomain 
in the x direction and N y = 32 interpolation points in the y direction. The table gives, 
in addition to the values of <x* X) > and Ax*, the asymptotic value x* s = 1— em(2/<5o). 
The average location of the transition layer is predicted very well by the asymptotics. 

Table 6 

Location of the transition layer at steady state; boundary data (UaJ- 



AS 


<«So > 




Ax* 


0.25 ■ 10"* 


0.40808 


0.40811 


0.27278 ■ 10"* 


0.50 ■ lO"* 


0.41277 


0.41241 


0.37231 • 10"* 


1.0 ■ 10"* 


0.42096 


0.42052 


0.92263 • 10 - * 


2.0 ■ 10~* 


0.43573 


0.43506 


1.4589 ■ 10~* 


3.0 ■ 10"* 


0.44854 


0.44782 


2.3924 • 10~* 



Again, the maximum deviation Ax* grows with AS, although not linearly as in the 
case of the step boundary data . 

Table 7 shows the impact of grid refinement (the number of collocation points per 
subdomain, N x , and the number of discretization points, N y ) on the computed value 
of x% Q {y) for the problem with boundary data (fll|), AS = 0.01. 

Table 7 

Effect of arid, refinement (N x , N y ) on <xj > (upper entries) and Ax* (lower entries); bound- 
ary data mi) with AS = 0.01. 



Ny 


N x = 19 


N x = 29 


N x =39 


N x =49 


N x =59 


8 


0.47848 
4.7396 10-* 


0.47653 
4.8375 10" 2 


0.47523 
4.8902 10" 2 


0.47557 
5.2615 10" 2 


0.47547 
5.2251 10" 2 


16 


0.47807 
4.7396 10~ 2 


0.47614 
5.8062 10" 2 


0.47610 
4.8849 10" 2 


0.47478 
4.9273 10" 2 


0.47516 
4.9518 10~ 2 


32 


0.47883 
4.7396 10-* 


0.47611 

4.8374 10" 2 


0.47578 
4.8864 10" 2 


0.47536 
4.9267 10" 2 


0.47481 
4.9516 10~ 2 


64 


0.47847 
4.7396 10-* 


0.47610 
4.8372 10" 2 


0.47561 
4.8867 10" 2 


0.47513 
4.9266 10" 2 


0.47524 
5.2280 10~ 2 



An obvious way to parallelize the algorithm ( f40| ) is to partition the interval (— n, ir) 
into subintervals of equal length 2tt/N v and map this partition onto a ring of proces- 
sors. Thus, one can achieve high speedups on a Paragon using nonblocking commu- 
nications. If each processor covers at least four mesh points in the y direction, only 
nearest- neighbor communication is needed. Table 8 gives ample evidence that the 
algorithm ( ]40|) is highly scalable; doubling the number of processors with the problem 
size results in almost identical CPU times. 

Table 8 

CPU time for 1,000 time steps on the Paragon XP/S as a function of the number of processors 
(P) and the size of the problem (measured by N y ); N x = 49. 



Ny 


P = 1 


P = 2 


P = 4 


P = 8 


P = 16 


32 


129.16 


65.88 


33.91 


17.43 




64 


252.85 


130.02 


65.85 


33.85 


17.45 


128 


519.05 


257.85 


129.68 


65.85 


33.92 


256 




515.01 


257.49 


129.77 


65.93 



Additional parallelism can be introduced by decomposing the domain in the x 
direction. However, our experience with a similar algorithm for combustion problems 
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indicates a potentially significant decrease (as much as 70%) in the efficiency of the 
algorithm ]35j| . 

In general, the algorithm ( fi"o| ) is very well adapted to the quasi one-dimensional 
structure of the transition layer. The algorithm predicts the location of the transition 
layer at steady state with a significant accuracy. The time step is of the same order 
of magnitude as for the one-dimensional analog ([))). 

3.4. Long-Time Integration. The algorithm ( |40| ) needs to be modified for 
long-time integration. We distinguish between the cases (3 = and (3 ^ 0. 

If (3 = 0, we use an algorithm similar to the one described in § 2.4. We start the 
integration of Eq. ( |2(i| ) at t = to, say. We identify the point x* with the location of 
the zero of the approximation U of u, averaged over y, at t = to and define the profile 
function uq as in Eq. (SF 



x — x 

(43) uq(x) = — tanh — , 

Then we integrate the nonlinear boundary- value problem 

— eAv + v t + Uf)V x + u'qV + 5o£~ 2 vv x = on (—1, 1) x (— tt, n) 

forward in time, subject to the boundary conditions 

«|«=-i = 6^e 2 (l + 6 - uo(-l)), «| B= i - ^V(-l - «q(1)), 

using the algorithm 

yn _ yn-l 

(44> eD 2 x V n + — + u D x V n + u' V n = eD 2 y V n ~ x ~ Sqe^V"- 1 D x V n - x , 

for n = 1,2,.... We define the approximation U of it, 

(45) U{x, y, t) = u (x) + 5o£~ 2 V{x, y, t), 

and integrate as long as the supremum of U{- , • ,t) —uq remains of the order of 6qe~ 2 . 
When this criterion is no longer met, at t = t± say, we suspend the integration, identify 
the point x* with the location of the center of the transition layer (averaged over y), 
and update the profile function uq. We repeat the procedure until the steady state is 
reached. 

The time step for the algorithm ( fy| ) is limited by the (explicit) term eDyV 71 ^ 1 , 
At < c(2n I Ny) 2 1 e , for some constant c < |. This limitation is not too severe, as e is 
very small and the variation of the solution in the y direction is exponentially small, 
so N y need not be large. 

If (3 ^ 0, the situation becomes more complicated. One can, of course, extend the 
algorithm (Q) trivially by incorporating the convective term in the right member, 

yn _ yn-l 
-eD 2 x V n + — + u D x V n + u' V n 



(46) = eD 2 yV n - x - Soe^V"- 1 DxV"- 1 - f3u Q D v V n - 1 - PSqe^V 71 - 1 DyV 71 ' 1 , 
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for n = 1, 2, . . . . Representative results obtained in this way for the boundary- value 
problem ( po| ) with f3 = 1 and e — 0.02 are given in Table 9. The boundary data are 
again piecewise constant, as in Table 5, but with So — 1.0 ■ 10~ 6 , 

( 1 + 1.0 ■ 10~ 6 - AS if - 7T < y < -lie, 
(47) u(-l,y)=< 1 + 1.0 ■ 10- 6 + A<5 if - \-n < y < ±tt, 

{ 1 + 1.0 ■ 10~ 6 - AS if±7r<y<7r. 

The algorithm (Q) was applied with N x = 39 collocation points per subdomain in 
the x direction and only N y — 16 interpolation points in the y direction. We observe 

Table 9 

Location of the transition layer at steady state; boundary data (U_J). 



A<5 




Ax* 


1.0 ■ 10~ b 


0.7099 


< 1.0 ■ 10" 4 


1.0 ■ lO -5 


0.7099 


< 1.0 ■ 10" 4 


1.0 • 10~ 4 


0.7101 


3.9 ■ 10~' 5 


0.5 ■ 10- J 


0.7122 


1.9 ■ 10-* 


1.0 • 10" a 


0.7162 


3.0 ■ 10-* 



that the average location of the transition layer does not change appreciably as long 
as AS is of the same order as the average perturbation So. As AS increases, the 
perturbation is no longer small compared with So, and the asymptotic results of § 3.1 
do not necessarily apply. Indeed, Ax* does not appear to vary linearly with AS, as 
was the case in Table 6. 



3 - 




-2 - 



-3|- ! ! ! ! ! t ! HI - 

-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 

x 

Fig. 2. Contour lines of the solution U at steady state. 

A graphical representation of the computed solution U and its deviation from the 
profile function, U — uq, for the case AS = 1.0 ■ 10~ 3 are given in Fig. 2 (contour lines 
of U at steady state) and Figs. 3 and 4 (perspective drawings of U and the difference 
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Fig. 4. The difference U — uo at steady state. 



U — u at steady state). The data for these figures were taken after 2,500,000 time 
steps (At = 0.4). 

The results are quite good, but for long-time integration one would do better by 
looking for v in terms of its Fourier coefficients Vk and adopting an implicit scheme 
in the y direction. Thus, the constraint on the time step becomes the same as in the 
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one-dimensional case. The algorithm would be of the following type: 

yn _ yn-l 

(48) eDlV? + k M k + u Q D x V£ + (i(3ku a + sk 2 )V£ + u' V n = -dos^W^ 1 , 

for n = 1, 2, . . . ; D y is the matrix of differentiation with respect to y in Fourier space, 
and Wk is some approximation to the kth Fourier coefficient of vv x + /3vv y . 

This algorithm parallelizes with respect to the Fourier modes. It has been applied 
to a problem involving a propagating combustion front in a moving fluid |35|, |3(| . 
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